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ABSTRACT 

The ongoing characterization of hot Jupiters has motivated a variety of circulation models of their 
atmospheres. Such models must be integrated starting from an assumed initial state, which is typically 
taken to be a wind-free, rest state. Here, we investigate the sensitivity of hot- Jupiter atmospheric 
circulation models to initial conditions. We consider two classes of models — shallow-water models, 
which have proven successful at illuminating the dynamical mechanisms at play on these planets, and 
full three-dimensional models similar to those being explored in the literature. Models are initialized 
with zonal jets, and we explore a variety of different initial jet profiles. We demonstrate that, in 
both classes of models, the final, equilibrated state is independent of initial condition. Otherwise 
identical models initialized with vastly different initial conditions all converge to the same statistical 
steady state. In some cases, the models exhibit modest time variability; this variability results in 
random fluctuations about the statistical steady state, but we emphasize that, even in these cases, 
the statistical steady state itself does not depend on initial conditions. Although the outcome of 
hot- Jupiter circulation models depend on details of the radiative forcing and frictional drag, aspects 
of which remain uncertain, we conclude that the specification of initial conditions is not a source of 
uncertainty, at least over the parameter range explored in most current models. 

Subject headings: hydrodynamics - methods: numerical - planets and satellites: atmospheres - planets 
and satellites: individual (HD 189733b, HD 209458b) 



1. INTRODUCTION 

Since the discovery of the fi rst exoplanet around a 
main-sequence star, 51 Pegasi b QMayor fc Queloz|1 995), 



tected. 



et al. 



2011) 



(|Wri£ 




Iau 


of 



giant planets orbiting within ~0 
their central stars. The physical regime of hot Jupiters 
differs substantially from those of solar-system giant 
planets; the presumed tidal locking leads to modest ro- 
tation rates and permanent day and nightsides, with in- 
cident stellar fluxes ~10 3 -10 6 times stronger than that 
received by giant planets in our solar system. These 
differences have motivated a flourishing research pro- 
gram focused on elucidating the atmospheric dynamics 
of these worlds. Current observations place meaningful 
constraints on planetary radii, atmospheric composition, 
albedo, and the three-dimensional temperature struc 
ture, including dayside temperature profiles and vari 
at ion of the temperature between day and night (e.g. 



Knutson et al.||20071 |Charbonneau et al.|[2 008; Knutson 
et al ||2008| |Cowan et al ||2007| |Harrington et al.| |2006, 



2007). 
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These observations have motivated a variety of three- 
dimensional (3D) circulation models of hot Jupiters 
dShowman fc Guillot| [20021 |Cooper fc Showman) [2005" 



2006 ; IShowman et ai.||2008U2009[ |Menou & Rauscher 



2009; Rausch eTTc Menou||2Q10| |2012| |Dobbs-Dixon ~& 
Lin||2008| |Dobbs-Dixon et al.||2010| |Lewis et al.||2010 
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Heng et al. 2011b a; Showman et al. 12012 ). These models 



obtain a circulation pattern generally comprising several 
broad atmospheric jets, including in most cases a strong 
eastward equatorial jet (superrotation) with speeds up 
to a few krns" 1 and westward zonal-mean flow at high 
latitude. An analytical explanation was provided by 



Showman & Polvani (2011), who showed that the super 
rotation results from the interaction of planetary-scale 
Rossby and Kelvin waves — themselves a response to the 
day-night thermal forcing — with the mean flow. When 
the radiative and advection time scales are comparable, 
the hottest region can be displaced eastward fro m the 
substellar point by te ns of degrees longitude (e.g., Show- 
man & Guillot 2002), as subse quently observed on HD 



189733b flKnutson et al.||2007[ ). Nevertheless, it is no- 



table that some 3D models produce qualitatively differ 
ent circulation patterns exhibi ting a range of flow behav- 
ior flThrastarson fc ChopOlO ). 

An important issue m modeling atmospheric circula- 
tion is the choice of initial conditions. Most of the models 
published to date integrate the equatio ns starting from a 
rest state co ntaining no winds, although Cooper & Show- 



man 



( 2005[ ) also performed an integration starting from 
an initial condition containing a broad westward jet and 
found that the initial con dition did not strongly affect the 
outcome. Nevertheless, Thrastarson & Cho (2010) re- 
cently performed a detailed investigation and found that, 
in their model, the initial conditions severely affect the 
final state. They explored initial conditions containing 
a broad eastward jet, broad westward jet, and three-jet 
pattern, all with speeds of 0.5-1 km s _1 , and compared 
this to models integrated from rest. These models — 
differing only in the initial condition — led to a wide range 
of final states whose pattern of flow streamlines and the 



2 



horizontal temperature distribution (including the lon- 
gitudinal o ffsets of any hot or c old sp ots) differed sig- 
nificantly. |Thrastarson fc Cho (2010) found that such 
extreme sensitivity to initial condition occurred even in 
models whose initial conditions differed only slightly. 

From a mathematical perspective, the initial condi- 
tions of course comprise an essential aspect of defining 
the overall mathematical problem, and in many dynam- 
ical systems, the initial conditions indeed play a crucial 
role in affecting the solution. However, atmospheres are 
forced-dissipative systems, and such forcing and dissipa- 
tion tend to drive the atmospheric circulation into a sta- 
tistical steady state retaining little if any memory of its 
initial condition. For this reason, many terrestrial gen- 
eral circulation model (GCMs) investigating the statisti- 
cal steady state of Earth's global atmospheric circulation 
are integrated from rest with no expectation that this 
choice adversely influences the results. 

A possible exception to this expected lack of sensi- 
tivity would be if the atmospheric circulation exhibits 
multiple, stable equilibria corresponding to radically dif- 
ferent circulation patterns for an identical set of forc- 
ing and boundary conditions. In this case, determining 
which of these multiple equilibria the atmosphere resides 
in requires knowledge of its history, potentially includ- 
ing its initial condition. An example in the context of 
terrestrial-planet climate is the existence of stable equi- 
libria corresponding to primarily ic e- free and glo bally 
glaciated ("Snowball Earth") states QBudyko 1969[ ). In 
the range of incident stellar fluxes where both equilibria 
exist, the actual state occupied by the climate depends 
on history (for speci fic examples of how this works, see 
Pierrehumbert 2010). In just such a way, the question 
ot initial-condition sensitivity for the atmospheric circu- 
lation is therefore perhaps best thought of as a question 
of whether the atmosphere exhibits multiple rather than 
only one stable equilibrium. To date, no claims have 
been made in the literature that the atmospheric circu- 
lation of hot Jupiters exhibit multiple, stable equilibria 
for a given set of forci ng conditions, and so a se nsitiv- 
ity to initial conditions ( |Thrastarson fc Cho|2010j ) is not 
expected a priori. Nevertheless, the issue is worthy of 
further study. 

Here, we present the results of a thorough exploration 
of the initial condition sensitivity in atmospheric circu- 
lation models of hot Jupiters. We investigate two classes 
of models — idealized shallow- water models in Section [2j 
as these have recently been used to identify dynamical 
mechanisms operating in th e atmospheres of hot Jupiters 
(Showman & Pojyani|2011 ) — and fully 3D models in Sec- 



tion [3J Section 3] concludes. 

2. SHALLOW WATER MODEL 

2.1. Model description 

Simplified models play an important role in under- 
standing atmospheric dynamical processes. The shallow- 
water model, in particular, is an idealized fluid sys- 
tem that has proven successful in illuminating many as- 
pects of the large-sc ale dynamics, both for Earth (e.g., 
Polvani et al. 1995) and giant planets (Dowling & In- 



4 This of course is not true for weather prediction models, which 
consist oTv ery-short-term integrations of typically a few da ys^ de" 



termining how some given observed circulation evolves with time. 




appropriate to synchronously 
Jupiters, the model produces a circulation with many 
similarities to those emerging in full 3D models of hot 
Jupiters. Here, we explore the sensitivity of this model 
to initial conditions. 

We adopt a two-layer model, with constant densities in 
each layer; the upper layer, of lesser density, represents 
the meteorologically active atmosphere, while the lower 
layer, of greater density, represents the deep interior. In 
the limit where the lower layer is quiescent and infinitely 
deep, this system reduces to the shallow- water equations 
for the flow in the upper layer: 



dv 

~dt 

dh 



■ gVh + /kxv = R- 



7"drag 



V • (vh) 



freq(A,0) - h 



7"rad 



(i) 



(2) 



where v(A, is horizontal velocity, h(\,(j),t) is the 
upper layer thickness, t is the time, g is the (reduced) 
gravity, / = 2Qsm(j) is the Coriolis parameter, k is 
the upward unit vector, Q is planetary rotation rate, 
d/dt = d/dt + v • V is the material derivative (including 
curvature terms in spherical geometry), and A and (j) are 
the longitude and latitude, respectively. 

Radiative heating and cooling are treated using a New- 
tonian cooling scheme, which relaxes the upper layer 
thickness toward a specified radiative-equilibrium thick- 
ness, /i e q(A, (j)) over a specified radiative timescale, r ra d. 
To represent the day-night heating pattern on a syn- 
chronously rotating hot Jupiter, h eq is chosen to be thick 
on the dayside and thin on the night side: 

(3) 



h eq = H + Ah eq cos A cos c 



where H is a constant mean thickness and Ah eq is the 
day-night contrast in radiative equilibrium thickness. 
The substellar point is at longitude 0° and latitude 0°. 
The equations include frictional drag, represented as a 
linear damping of winds with a specified drag timescale, 
which could represent the potential effec ts of magnetohy- 
drodynamic friction (|Perna et al.J|201Q[ ), vertical turbu- 
lent mixing (|Li & Goodman 2010), or momentum trans- 



port by breaking gravity waves TWatkins & Cho||2010b 



The term R in Eq. pj represents momentum transport 
between the layers ana is 



R(A,&t) 



0, 



Q>0; 
Q < 



(4) 



Air moving into the upper layer (Q > 0) affects the upper 
layer's specific angular momentu m, but air moving out 
of the upper layer does not. See Showm an fc Po lvani 
(2011) for further discussion and interpretation of the 



equations. Parameters were chosen to be representative 
of hot Jupiters, including g = 20 msec -2 , H = 200 km, 
ft = 3.2 x 10 -5 sec -1 and a = 8.2 x 10 7 m, similar to the 
values on HD 189733b. Note that the specific choices of 
these parameters are not essential to the result; similar 
results would obtain were other parameter values used 
instead. 
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Fig. 1. — Initial jet profiles adopted in our shallow-water models, 
corresponding to a rest state (black line), eastward jet (blue line), 
westward jet (green line) and three-jet pattern (red line). 



Thras- 



Our initial conditio ns are motivated by those ofl 
tarson & Cho (2010) and are shown in Figure |T] The 
model is initialized with a zonally symmetric (i.e., 
longitude-independent) zonal flow given by 



^initial ((/>) = U exp 



2a 2 



(5) 



where a is the jet half-width and U is the jet speed. 
We present models whose initial zonal winds comprise 
eastward equatorial jets (U = lkms -1 and a = 7r/12, 
blue curve in Figure m), westward equatorial jets (U = 
— lkms -1 and a = 7T/12, green curve in Figure [l]), and a 
three-jet profile (red curve), as well as models integrated 
from rest (black curve). The initial meridional velocity 
is set to zero and the height field is specified to be in 
gradient-wind balance with the zonal-wind field. Note 
that this is only a subset of the models explored; other 
initial conditions were tried as well and yielded the same 
results. 

We integrate Eqs. Q-pl) using t he Spectral Transform 
Shallow Water Model (STSWM) (jHack fc Jakob |[l992| , 
which solves the equations using pseudospectral methods 
in vorticity-divergence form. We adopt a spectral trun- 
cation of T170, corresponding to a resolution of 0.7° in 
longitude and latitude (i.e., a global grid of 512 x 256 in 
longitude and latitude). The code adopts the leapfrog 
time stepping method, using an Asselin filter to suppress 
the computational mode. A V 6 hyperviscosity is applied 
to each of the dynamical variables to maintain numeri- 
cal stability. All models are equilibrated to a statistical 
steady state. 

2.2. Results 

We explored a wide range of models with differing val- 
ues of r ra d, Tdrag, Ah eq _/ H , and initial condition. A small 
subset of these models, which are illustrated in subse- 
quent figures, are shown in Table [T] 

We find, over a wide range of parameters, that the 
final, equilibrated state is independent of initial condi- 
tion. This is illustrated in Figures [2] and [3} At suffi- 
ciently low amplitude — that is, at sufficiently low values 
of Ah eq /H — the equilibrated solutions are temporally 
steady, whereas time variability sets in beyond a critical 



amplitude ( [Showman fc Polvani| [20101 12011[ ). Figure [5] 
shows the steady-state geopotential, gh\ and zonal- mean 
zonal winds at an integration time of 100 day^] in low- 
amplitude models with Ah eq /H = 0.01, r ra d = 1 day, 
and Tdrag = 1 day. The four models in Figure |2| are inte- 
grated from the four initial conditions shown inFigure[T] 
corresponding to a rest state (top row), eastward equa- 
torial jet (second row), westward equatorial jet (third 
row), and three-jet pattern (fourth row). As Figure [2] 
demonstrates, all aspects of the equilibrated, steady- 
state flow field — including the spatial pattern of the equi- 
librated geopotential and the zonal-mean zonal winds — 
are essentially identical regardless of the initial condition 
used. This final state consists of standing, planetary- 
scale Rossby and Kelvin waves; two anticyclones straddle 
the equator on the dayside and two cyclones s traddle it 
on the nightside ( see | Showman fc Polvani|2011| ). Because 
of the low forcing amplitude, the equilibrated zonal- mean 
zonal wind is weak — corresponding to an equatorial su- 
perrotating flow with a speed of only 0.01ms -1 . Note 
that all models equilibrate to this identical final jet pro- 
file despite the fact that the speed of the initial jet, 
±1 km s -1 , exceeds that of the equilibrated jet by a factor 
oflO 5 . 

Figure [3] illustrates an example of a high- amplitude 
model, where Ah eq /H = 0.5, r ra d = 0.1 day, and 
Tdrag — » oo (meaning there is no ex plicit large-scale drag 
in the upper layer; as described in |Showman fc Polvani 



( |2011[ ), such a model still equilibrates because of inter- 
actions with the quiescent lower layer). Again, the four 
models in Figure [3] are integrated from the four initial 
conditions shown in Figure [I] corresponding to a rest 
state (top row), eastward equatorial jet (second row), 
westward equatorial jet (third row), and three-jet pattern 
(fourth row). All of the models equilibrate to the same fi- 
nal state, with significant day-night differences in geopo- 
tential and an overall pattern of eastward-equatorward 
phase tilts, particularly on the dayside, which is the re- 
sult of the standing, planetary-scale Rossby and Kelvin 
waves. Because of the large forcing amplitude, short r ra d, 
and absence of large-scale drag, the zonal-mean zonal 
wind equilibrates to fast speeds of 1 kms -1 in the core of 
the superrotating equatorial jet that emerges (Figure [3]). 
Again, we emphasize that the speed and amplitude of 
this equilibrated jet is totally independent of whether the 
initial condition contained an eastward jet, a westward 
jet, multiple jets, or no jets at all. 

Figure [4] further quantifies the similarity between 
the final states of otherwise identical models initial- 
ized with differing initial conditions. The left pan- 
els show the differences in the geopotential, at a 
given time, in the equilibrated states of two mod- 
els integrated with identical forcing parameters but 
differing initial conditions, i.e., [gh mo d e \ a (\, — 
b(A,</>, ti)]/gh mode \ a (\,(j),ti)), where "model a" 
and "model b" are the two models being compared, and 
t\ is some late time after the runs are equilibrated. The 
right panels overplot the zonal-mean zonal wind for these 
two models in red and green. The top row represents 
the differences between two low-amplitude models (SW2 
and SW3) while the bottom row shows the differences 
between two high- amplitude cases (SW5 and SW8). 



5 In this paper, 1 day is defined as 86400 sec. 
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TABLE 1 

The properties of some shallow- water runs 



Name 


Initial condition 


Aheq/H 


Trad (days) 


Tdrag (days) 


Time (days) 


SW1 


Rest state 


0.01 


1 


1 


100 


SW2 


Eastward jet 


0.01 


1 


1 


100 


SW3 


Westward jet 


0.01 


1 


1 


100 


SW4 


Combined jet 


0.01 


1 


1 


100 


SW5 


Rest state 


0.5 


0.1 


oo 


1000 


SW6 


Eastward jet 


0.5 


0.1 


oo 


1000 


SW7 


Westward jet 


0.5 


0.1 


oo 


1000 


SW8 


Combined jet 


0.5 


0.1 


oo 


1000 
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Fig. 2. — Geopotential (gh) over the globe (left) and zonal- mean 
zonal wind (right) in four low-amplitude shallow-water models ini- 
tialized with differing initial conditions corresponding to the four 
initial conditions in Figure ^ a rest state (top row), eastward 
equatorial jet (second row), westward equatorial jet (third row), 
and three-jet pattern (fourth row). Units of colorbar are m 2 s -2 . 
These are models SW1, SW2, SW3, and SW4 from top to bot- 
tom, respectively. These are snapshots, shown at 100 days, after 
the models have reached equilibrium. Despite the different initial 
conditions, the models all converge to the identical steady state. 

When the forcing amplitude is low and the solutions 
are steady, the final solutions are identical, in a point-to- 
point sense, to a precision of literally ~10 -12 -10 -13 (top 



2.5 



3.5 



4.5 



x10 D 



Fig. 3. — Geopotential (gh) over the globe (left) and zonal-mean 
zonal wind (right) in four high-amplitude shallow-water models ini- 
tialized with differing initial conditions corresponding to the four 
initial conditions in Figure ^ a rest state (top row), eastward 
equatorial jet (second row), westward equatorial jet (third row), 
and three-jet pattern (fourth row). Units of colorbar are m 2 s -2 . 
These are models SW5, SW6, SW7, and SW8 from top to bot- 
tom, respectively. These are snapshots, shown at 1000 days, after 
the models have reached equilibrium. Despite the different initial 
conditions, the models all converge to the same statistical steady 
state. 



row of Figure [4|. Fractional differences are a few xlO 13 
over most of the globe but rise to a fewxlO -12 in a 
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few localized regions (particularly near the poles) . These 
miniscule differences are numerical, resulting from a com- 
bination of roundoff and discretization error, and indi- 
cate that, for all practical purposes, the solutions of these 
different models are truly identical despite the vastly dif- 
ferent initial conditions. The equilibrated zonal- mean 
zonal wind profiles are likewise so similar that the red 
curve is precisely covered by the overlying green curve 
(top right panel of Figure [4]). 




500 1000 

U [msec -1 ] 

Fig. 4. — Differences between the models in Figures [2] and [3] 
Left: Fractional differences in geopotential (gh) between two oth- 
erwise identical models initialized with different initial conditions. 
Top left plots (g/isw2 _ ^sw3)/fi , ^SW2 J whereas bottom left plots 
(ghsW5 ~ 5^SW8)/^SW5i both at a specific instant in time (100 
days in the top row, 1000 days in the bottom row). Here, ghx 
refers to the geopotential of model X. Right: Zonal-mean zonal 
wind from the two models shown in the corresponding left panels. 
Top right shows SW2 (green curve) and SW3 (red curve) at 100 
days while bottom right shows SW5 (green curve) and SW8 (red 
curve) at 1000 days. The figure demonstrates that low-amplitude 
cases (SW2 and SW3) converge to the identical steady state, with 
fractional differences of ~10 -12 or less; high-amplitude cases, how- 
ever, converge to states that, at any given instant, can differ at spe- 
cific points by up to ~2%, although the zonal-mean zonal winds 
are still nearly identical. 

Although Figure [3] demonstrates that high-amplitude 
models likewise equilibrate to nearly identical states, tiny 
differences become apparent when one compares in more 
detail. The bottom row of Figure H] quantifies these dif- 
ferences. In this case, the fractionalpoint-to-point differ- 
ences at a given time between otherwise identical models 
initialized with differing initial conditions (in this case 
SW5 and SW8) reaches ~1% in specific regions, particu- 
larly on the night side. One can ask whether this is a true 
difference in the statistical steady states between these 
models or whether it rather results from time variability 
that might induce a slight randomness around a single 
statistical steady state. To address this question, Fig- 
ure [5] shows the differences between two different snap- 
shots at different times within a given model integration. 
In the bottom row, these differences are likewise seen 
to reach ~1%, with a spatial pattern extremely simi- 
lar to that seen in the bottom row of Figure |4j This 



comparison indicates that the model-model differences 
shown in Figure [4] are the simple result of time variability 
and not the result of any fundamental sensitivity of the 
statistical steady state to initial condition. Indeed, the 
differences between the time- averages of the same two 
models are much smaller than the differences between 
their instantaneous snapshots displayed in Figure [4] (not 
shown), confirming that the statistical steady states are 
essentially identical for these models. We also note that, 
even at a given time, the two models have nearly identi- 
cal zonal- mean zonal wind profiles (Figure |4j lower right 
panel), indicating that even the instantaneous point-to- 
point variability has little effect on the zonal-mean state. 

We emphasize that the results shown are not spe- 
cific to the particular parameters illustrated in the fig- 
ures but rather are general. The wide variety of models 
we have performed all confirm the essential point made 
here, namely, the insensitivity of this forced shallow- 
water model to initial conditions. 




500 1000 

longitude u [msec" 1 ] 

Fig. 5. — Differences between a given high-amplitude shallow- 
water model, SW6, at two different points in time (900 and 
1000 days), after the models have reached statistical steady 
state. Left panel shows the fractional geopotential differences, i.e., 
(s^900days - 5 f ^iooodays)/5 f ^900days, whereas right panel shows the 
zonal-mean zonal wind at 900 days (green curve) and 1000 days 
(red curve). These temporal differences within a given model are 
similar to the model-model differences shown in the bottom row 
of Figure ^] demonstrating that those differences result from time 
variability rather than sensitivity of the statistical steady state to 
initial conditions. 



3. THREE-DIMENSIONAL MODEL 

3.1. Model description 

We now consider 3D models of the atmospheric cir- 
culation. As in most previous investigations, we adopt 
the primitive equations. We solve the equations us- 
ing the MITgcm, which is a state -of-the-art circulation 
model flAdcroft et aLl|2QQ4[ ) that |Showman et al.| ( [2009 ) 
adapted tor application to hot Jupiters. The horizontal 
momentum, vertical momentum, continuity, and thermo- 
dynamic energy equations are, using pressure as a verti- 
cal coordinate, 



- V$ -/kxv + P v 

d§ _ 1 

dp p 

duo 
dp 

dT _ q uo 
dt c p pc p 



V- v 







(6) 
(7) 
(8) 
(9) 
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where v is the horizontal velocity on constant-pressure 
surfaces, uj = dp /dt is the vertical velocity in pressure co- 
ordinates, <I> is the gravitational potential on constant- 
pressure surfaces, / = 2Q sin cj> is the Coriolis parame- 
ter, Vt is the planetary rotation rate, k is the local ver- 
tical unit vector, q is the thermodynamic heating rate 
(Wkg -1 ), and T, p, c p are the temperature, density, 
and specific heat at constant pressure. V is the hori- 
zontal gradient evaluated on constant-pressure surfaces, 
and d/dt = d/dt + v • V + cod /dp is the material deriva- 
tive (including curvature terms in spherical geometry). 
The term V v is a velocity damping term, including a 
Shapiro filter to maintain numerical stablity (which has 
only a small effect on the large-scale flow), and option- 
ally, an explicit large-scale frictional drag term (see be- 
low). Eq pi is actually solved in an alternate form, 



d0 
dt 



T c n 



(10) 



where = T(p/po) K is the potential temperature (a mea- 
sure of entropy) , k is the ratio of gas constant to specific 
heat at constant pressure, and po is a reference pressure 
(note that the dynamics are independent of the choice of 
Po). The dependent variables v, oj, <£, p, and T are 
functions of longitude A, latitude 0, pressure p and time 
t. 
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Fig. 6. — Forcing profiles used in the Newtonian heating/cooling 
scheme. Top: Radiative-equilibrium temperature profiles. Dashed 
curves show the radiative equilibrium temperature profile for the 
nightside (T n i g h t eq , left curve) and substellar point (right curve). 
The solid curve shows Ti ro (p). The dotted curve shows the 
substellar-nightside radiative-equilibrium temperature difference, 
ATeq(p), assuming the temperature difference at the top, AT eq ,top, 
is 1000 K. Bottom: Radiative time constant, r ra d, versus pressure. 



Showman et al.J (2009]) coupled the MITgcm to the 
multi-stream radiative transfer model of Marl ey fc| 



heating rates when the atmospheric composition and 
opacities are specified. In the present context, however, 
our goal is to characterize the sensitivity to initial condi- 
tions in the clearest possible context, and so rather than 
using this coupled model, we specify the radiative heat- 
ing/cooling using a Newtonian cooling scheme, which 
relaxes the temperature toward a specified radiative- 
equilibrium temperature over a specified time constant: 



T(\,</),p,t) -r eq (A,<ft,p) 

Tradb) 



(11) 



The Newtonian co oling scheme has been widely used in 
exoplanet studies dShowman fc G uillot 20 02] |Cooper fc 
Showman 2005; Showman et al. 2008; Menou & Rauscher 



2009; Rausch erir Menou|2010l|Perna et al.|2010||Thras- 
|tarson fc Cho|2010| |Heng et al.||2011bp . 

The radiative-equilibrium temperature, T eq (A,0,p), is 
defined as 



T eq (A,<M 



flight, eq(p) 
2~night,eqCp) J 



AT, 



where T nightjeq (p) + AT eq is 



e q{p) cos A cos </>, dayside; 

nightside 
(12) 

the radiative-equilibrium 



McKay ( |1999[ ), which allows for accurate calculation of 



temperature at the substellar point and AT eq (p) is the 
difference in radiative-equilibrium temperature between 
the substellar point and the nightside. As written this 
expression takes the substellar point to be at longitude 
and latitude of 0°,0°. To specify the nightside pro- 
file, we further define T n i g ht, e q(p) = ^iro — A ^ eq , where 
Tixoip) is the one-dim ensional radiativ e-equilibrium tem- 
perature profile from Iro et al.] ( |2005| ). These definitions 
then imply that the substellar radiative-equilibrium tem- 
perature profile is T\ TO + ^^ SL . Since radiative heating 
is strong at the top and weak at the bottom, it is im- 
portant that AT eq (p) decreases with increasing pressure. 
For computational simplicity we specify AT eq (p) as a 
piecewise-continuous analytical function that is a con- 
stant, AT eq5t0 p, at pressures less than p e q,top; is zero at 
pressures exceeding p e q,bot5 and varies linearly with log- 
pressure in between. For the models described in this 
paper, we take AT eq5top = 1000 K, p eq ,to P = 10" 3 bar, 
and Peq,bot = 10 bar; note that our key result — namely, 
insensitivity to initial conditions — does not depend on 
these precise values. The nightside and substellar radia- 
tive equilibrium profiles, as well as T\ TO and AT eq (p), are 
shown in Figure [6j 

Likewise, the radiative time constant is expected to be 
a strong function o f pressure, being short at the top and 
long at the bottom fllro et al.|2005{ Showman et al.|2008[ ). 
For computational simplicity, we here assume that r ra d 
is a function of pressure alone, and we again define a 
piecewise-continuous analytic function that allows such 
a downward-increasing behavior: we take r ra d(p) to be 
a constant, r ra d,top 7 at p < Prad,top5 another constant, 
Trad,bot, at p > Prad,bot; and we assume that logr rad 
varies linearly with logp in between. In this paper, we 
take Prad,top = 10 _2 bar and p ra d,bot = 10 bar. We ex- 
plore several values for r ra d,top and r ra d,bot in different 
models, with r ra d,top generally chosen to be short and 
7"rad,bot chosen to be long. Again, our key results are not 
dependent on the precise values. The profile of r ra d for 
one such model is shown in Figure [6] 
In our models, we also include a simple frictional drag 
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Fig. 7. — Initial jet profiles adopted in our 3D models. Plotted is zonal-mean zonal wind of the initial condition, in ms — 1 , versus latitude 
and pressure. Initial conditions correspond to a rest state (upper left), eastward decaying jet (upper right), westward decaying jet (lower 
left), and eastward barotropic jet (lower right). We also tried westward barotropic jets in some models (not shown). 



scheme near the bottom of the domain. This might 
crudely represent the effects of "magnetic drag" asso 
ciated with the partial ioniz ation expected at temper- 



ature s exceeding ~1500K (Perna et al. 2010 Menou 



small values at the bottom, as expected on real hot 
Jupiters; (2) to produce equilibrated circulation patterns 
qualitatively resembling those from mo dels that couple 
the dynamics to radiative transfer (i.e., | Showman et al. 



A bar (see Figure |6|. From a more practical perspec- 
tive, such frictional drag also forces the flow to equi- 
librate in a reasonable integration time. When study- 
ing sensitivity to initial conditions, it is particularly im- 
portant to ensure that the models have reached equilib- 
rium, and this is aided by including such a drag scheme. 
The drag is introduced on the righthand side of Equa- 
tion (|6| and takes the form — k v v, where k v (p) is a 
pressure-dependent drag coefficient. The drag coeffi- 
cient is zero at pressures less than £>drag,top and equal 

tO k F (p - Pdrag,top)/(Pbot -Pdrag,top) at p > Pdrag,top, 

where Pdrag,top is the lowest pressure of the region expe- 
riencing drag, pbot is the mean pressure at the bottom 
of the domain, and kp is a constant (t his formulation 
of dr ag is extremely similar to that of |Held fe Suarez] 
|1994| ). This formulation implies that the drag coefficient 
increases linearly with pressure from zero at Pdrag,top to 
kp at the bottom of the domain. Motivated by expecta- 
tions that magnetic drag is most important only at tem- 
peratures exceeding ~1500K, we take Pdrag,top = lbar 
in most models, although we also explore values of 10 bar 
to determine the sensitivity to drag scheme. The qualita- 
tive structure of the equilibrated dynamical state is not 
strongly sensitive to kp] here, we explore values of 0.1 
and 0.01 day -1 , implying characteristic drag timescales 
of 10 and 100 days near the bottom. 

Overall, our choices of AT eq (p), r ra d(p) and drag 
scheme described above are motivated by three over- 
arching goals: (1) to ensure that the radiative heat- 
ing/cooling rates (expressed in Ks -1 ) are large at the 
top but decrease rapidly with increasing pressure to very 



2012[ ), which occur in our model at pressures exceeding 2009; Heng et al. 201 1a] [Rauscher & Menou||2012| |Perna 

id 



et al. 2012| ), and (3) to ensure that the models equili- 
brate m finite time, as necessary to test sensitivity to 
initial conditions and to survey the parameter space. 
The first and second criteria generally lead to choices of 
^rad, top = 10 4 or 10 5 s and r ra d,bot ^ 10 6 s, while the third 
criterion suggests T ra d DO t 

< 10 8 s and 

Pdrag,top ^$ 10 bars. 

We no te that our fomulat i on di ffers significantly from 
that of iThrastarson fc Cho| ( |2010| ) , who choose AT eq and 
r ra d to be independent ol pressure; although this assump- 
tion has the advantage of simplicity, it fails to satisfy 
criteria (1) and (2). 

Following jThrastarson fc Cho| ( |2010| ) , as well as our 
shallow- water models from Section [2j we initialize most 
of our 3D models with a zonally symmetric zonal flow 
whose latitude dependence is given by Equation ([5| 
with a = 7r/9 and U = (corresponding to a rest 
state), lkms -1 (corresponding to an eastward equa- 
torial jet), or —lkms -1 (corresponding to a westward 
equatorial jet). In some cases, we assume this ini- 
tial jet to be independent of pressure, while in oth- 
ers, we allow the jet to decay with pressure by multi- 
plying the right side of Equation ^ by the function 
1.6+arctan[— (logp— logptop)/logPtopJ/l-6, which causes 
the jets to decay from a peak speed of ~1.6kms -1 at 
the top to ~400ms -1 at the bottom of the domain. Fig- 
ure [7] shows several of these initial conditions, laid out in 
a format that we will repeat, for easy comparison, when 
presenting results. 

We adopt planetary parameters appropriate to 
hot Jupiters, including specific heat c p = 1.3 x 
10 4 Jkg -1 K -1 , specific gas constant 3700 Jkg -1 K -1 , 
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TABLE 2 
The properties of some 3D runs 



Name 


Initial condition 


7"rad,top ( s ) 


frad,bot (s) 


kF (day" 1 ) 


Pdra S ,to P (bar)) 


GCM1 


Rest state 


10 4 


10 6 


0.1 


1 


GCM2 


Eastward decaying jet 


10 4 


10 6 


0.1 


1 


GCM3 


Westward decaying jet 


10 4 


10 6 


0.1 


1 


GCM4 


Eastward barotropic jet 


10 4 


10 6 


0.1 


1 


GCM5 


Rest state 


10 5 


10 6 


0.1 


1 


GCM6 


Eastward decaying jet 


10 5 


10 6 


0.1 


1 


GCM7 


Westward decaying jet 


10 5 


10 6 


0.1 


1 


GCM8 


Eastward barotropic jet 


10 5 


10 6 


0.1 


1 



and a rotation period, gravity, and planetary radius of 
3.024 x 10 5 s, 9.36 ms- 2 , and 9.437 x 10 7 m, respectively. 
These values are appropriate to HD 209458b, although 
we emphasize that our results are not sensitive to the pre- 
cise values, and similar behavior would occur had choices 
appropriate to other typical hot Jupiters been made in- 
stead. 

The MITgcm solves the equations using a finit e-volume 
discretization o n staggered Arakawa C grid ( [Arakawa] 
Lamb | |1977| ). Rather than the standard longi- 
tude/latitude coordinate system, we so lve the equations 
using the cubed-sphere grid following Showman et al 



(2009). The horizontal resolutions is C32 in most mod- 



els, implying that each of the six "cube faces" has a reso- 
lution of 32 x 32 finite- volume elements, which is roughly 
equal to a global resolution of 128 x 64 in longitude and 
latitude. However, to ensure that our results are numer- 
ically converged and do not depend on these numerical 
details, we also performed some models at a resolution of 
C64 (i.e., 64 x 64 cells on each cube face, corresponding 
to a global resolution of approximately 256 x 128) and a 
resolution of C128 (i.e., 128 x 128 on each cube face, 
corresponding to a global resolution of approximately 
512 x 256). The upper boundary is zero pressure and 
the bottom boundary is an impermeable surface. We 
adopt Nl = 40 levels in the vertical; the bottom Nl — 1 
levels are evenly spaced in log-pressure between 200 bars 
at the bottom and 0.2 mbar at the top; the top layer 
extends from a pressure of 0.2 mbar to zero. 

3.2. Results 

As before, we explored a variety of models, with dif- 
fering values of T ra d ? top? r rad,bot? 

kF, 

Pdrag,top an d initial 
condition. A small subset of these models, which are 
illustrated in subsequent figures, are shown in Table [2j 

In agreement with our shallow- water results, we find 
that the final, equilibrated state of our 3D models are 
independent of the initial condition. This is illustrated 
for a particular choice of forcing parameters (correspond- 
ing to models GCM1 to GCM4) in Figures M and M For 
easy comparison, the four panels in each of these figures 
adopt the initial conditions of the corresponding panels 
of Figure [7] — a rest state (top left panel), eastward decay- 
ing jet (top right), westward decaying jet (bottom left), 
and eastward barotropkj^] jet (bottom right). Figure [8] 
shows the zonal-mean zonal wind, while Figure [9] shows 
the temperature and two-dimensional velocity pattern at 
a pressure of 30 mbar, for these four models after equilib- 
rium has been reached. As can be seen, all four models 

6 In this context, "barotropic" means that the horizontal wind 
speed is independent of pressure. 



exhibit extremely similar patterns of zonal wind, tem- 
perature pattern, and two-dimensional velocity structure 
despite the differing initial conditions. The momentum 
fluxes caused by the day-night thermal forcing drive a 
broad equatorial jet whose peak speeds exceed 3kms _1 
(Figure pi). The day-night temperature differences ex- 
ceed ~8lJ0K at the top of the model and are 500 K at 
30 mbar (Figure [9]). At low pressure, the short radiative 
time constant — 10 4 s in these models — leads to little lon- 
gitudinal offset of the dayside hot region, although by 30 
mbar the hot spot is displaced to the east of the substel- 
lar point by ~20° longitude. Figure [l0| shows the total 
kinetic energy over time, integrated over the entire do- 
main, for these four models; the initial kinetic energies 
differ because of the differing initial jet profiles, but the 
models all converge to the same kinetic energy over time. 
Clearly, the models have lost memory of the initial con- 
dition and have all converged to the identical statistical 
steady state. 

The lack of sensitivity to initial conditions is demon- 
strated with another example in Figures 11 and[l2j which 
show models with a longer radiative time constant of 
10 5 s in the upper part of the domain. As before, the 
figures depict the zonal- mean zonal wind (Figure 11) 
and the temperature and velocity patterns at 30 mDar 
(Figure pk for four models— GCM5, GCM6, GCM7, 
and GCM8 — initialized with the jet profiles shown in 
Figure [7| Except for the initial conditions, everything 
about the models are identical. Despite the vastly dif- 
fering initial conditions, the models again all converge 
to the same final state. The equilibrated state exhibits 
a fast (3kms _1 ) equatorial jet, with day-night tempera- 
ture differences that are hundreds of K at low pressure. 
Interestingly, because of the longer radiative time con- 
stant at the top, the day-night temperature difference in 
GCM5-GCM8 is smaller than in models GCM1-GCM4, 
and there is a larger eastward displacement of the dayside 
hot region relative to the substellar point. Thus, while 
it is clear that the response of our hot- Jupiter models 
do depend on forcing parameters, they do not depend on 
initial conditions. 

If one compares the panels in Figures [TT] and [12] care- 
fully, very slight differences become apparent; the peak 
speed at the core of the equatorial jet, for example, are 
not quite identical between the four panels. Likewise 
there are slight, second-order fluctuations in the pattern 
of the velocity vectors in Figure [9j particularly on the 
nightside within ~40° latitude of the equator. These dif- 
ferences are the result of time- variability, which induces 
a slight randomness about the statistical steady state, 
rather than any fundamental sensitivity to initial con- 
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Fig. 8. — Equilibrated zonal-mean zonal wind in four otherwise identical models with differing initial conditions. These four models 
adopt the initial conditions corresponding to the equivalent panels of Figure [7] Upper left, upper right, lower left, and lower right panels 
show models GCM1, GCM2, GCM3, and GCM4, initialized from a rest state, an eastward decaying jet, a westward decaying jet, and an 
eastward barotropic jet, respectively. The zonal averages shown here are time averaged from 700 to 1000 days. All four models adopt 



T r ad,top = 10 4 s, T rad)bot = 10 6 s, k F = 0.1 day \ p d 



rag, top 



1 bar, and AT e 



eq,top : 



1000 K. 
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Fig. 9. — Temperature (color scale, in K) and winds (arrows) for the same four models shown in Figure [8] Panels show models GCM1 
(upper left), GCM2 (upper right), GCM3 (lower left), and GCM4 (lower right), initialized with a rest state, eastward decaying jet, westward 
decaying jet, and eastward barotropic jet, respectively. These are snapshots at 1042 days. The final states are extremely similar. 
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Fig. 10. — Total kinetic energy, integrated over the model do- 
main, for the four models GCM1 (black), GCM2 (green), GCM3 
(red) and GCM4 (blue) which are initialized with different initial 
conditions but are otherwise the same. The initial kinetic ener- 
gies differ greatly due to the differing initial conditions, but the 
model all converge to the same final state, independent of initial 
condition. 



ditions. This is demonstrated in Figure [T3j where the 
zonal-mean zonal wind from model GCM8 is shown at 
four different times after the model has reached statisti- 
cal steady state. The figure shows that the equatorial jet 
fluctuates slightly in time; the amplitude of these fluctu- 
ations is comparable to the inter-model differences seen 
in Figure [TT] A temporal average of the temperature 
or velocity patterns removes these random fluctuations 
and yields a pattern that is essentially identical between 
models with the same forcing parameters but differing 
initial conditions. 

We also performed models with weaker damping and 
drag in the deep atmosphere, for example models with 
Trad,bot of 10 7 s, kF of 0.01 day -1 , and/or p d rag,to P of 
lObar (not shown). We also performed some models 



ureQbut multiplied by —1) were explored. All of these 
models converged to a statistical steady state that is es- 
sentially identical, showing that memory of the initial 
condition has been lost. For a given set of forcing pa- 
rameters, the overall pattern of temperature and winds 
are also extremely similar between our C32, C64, and 
CI 28 models, suggesting that numerical convergence has 
been essentially reached even by C32. 

4. DISCUSSION AND CONCLUSIONS 

We explored the sensitivity to initial conditions of 
three-dimensional models of synchronously rotating hot 
Jupiters with day-night thermal forcing. The thermal 
forcing was chosen to be strong at the top and weak at 
the bottom, as must occur on real hot Jupiters. Models 
were integrated from rest and from various eastward and 
westward jet profiles with speeds up to ~1.5kms _1 . In 
all models explored, we found that the statistical steady 
states are independent of initial conditions; for a given 
set of forcing parameters, models converged to the same 
statistical steady state regardless of the initial condition 
employed. When the thermal forcing is strong, the cir- 
culation in the equilibrated state exhibits modest time 
variability that induces small-amplitude random fluctu- 
ations in any given realization. The statistical steady 
state itself, including not only the time-mean wind and 
temperature but the overall amplitude of these fluctua- 
tions, are independent of initial conditions in all of our 
models. 

As described in the Section [I] the issue of initial- 
condition sensitivity is perhaps best thought of in terms 
of whether the atmospheric circulation exhibits a single, 
rather than multiple, stable equilibria. Taken at face 
value, our models suggest empirically that, for any given 
set of forcing and damping parameters, there exists only 
one stable equilibrium — at least for the range of forcing, 
damping, and planetary parameters explored here. We 
have intentionally chosen forcing and planetary parame- 
ters similar to those appropriate to typical hot Jupiters, 
including HD 189733b and HP 209458b, as explored by 
a number of authors ( Showm an et al.| 2009 ; H eng et ah ] 
2011b|a| iPerna et al. 120121 IRauscher & Menou 201 



than 1000 K as for most of the models in this paper. 
Likewise, we also tried some qualitatively different ini- 
tial conditions, such as models containing nonzero zonal 
flow over only a specified subset of longitudes. In each 
case, for a given set of forcing/damping parameters, the 
models equilibrate to a statistical steady state with no 
sensitivity to initial conditions. 

The lack of sensitivity to initial conditions in our mod- 
els is not an artifact of our model resolution or numer- 
ical damping (i.e., the Shapiro filter) but rather is a 
fundamental property of the system behavior over the 
range explored. We also integrated models with reso- 
lutions of C64 and CI 28, approximately equivalent to 
global resolutions of 256 x 128 and 512 x 256, respec- 
tively, in longitude and latitude. To the best of our 
knowledge, our C128 models, in particular, are higher 
resolution than any published three-dimensional models 
of hot Jupiters that include day-night thermal forcing 
to date. Initial conditions corresponding to rest states, 
eastward barotropic jets, and westward barotropic jets 
(identical to that shown in the lower right panel of Fig- 



with weaker day-night forcing, e.g., AT eq = 100 K rather Therefore, we expect our fundamental result — ^e lack 

of sensitivity to initial conditions — to apply generally to 
the regimes explored in those papers. 

In this context, it is intere sting that our findings diffe r 
so drastically from those of |Thrastarson fc Choj pOlO). 
Their model differs from ours not only in its sensitivity 
to initial condition but also in that it seems to exhibit 
much greater temporal variability at essentially global 
scales. The reasons for these divergent behaviors are 
unclear but may result in part from the different ther- 
mal forcing in the two stu dies. In most of their models, 
Thrast arson & Cho (2010) adopt profiles of T eq and r ra d 
in their INewtoman-coolmg scheme that are independent 
of pressure, which means that their models exhibit large 
entropy variations against the impermeable (rigid) lower 
boundary. As is well known, the existence of horizon- 
tal entropy variations against an impermeable surface 
tend to make a flow much more baroclinically unstable 
(see, e.g., Vallis 2006, Chapter 6). Such instabilities can 
lead to significant time variability, particularly when the 
Rossby deformation radius is global in scale, and this 
may contribute to the large degree of temporal variabil- 
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Fig. 11. — Zonal-mean zonal wind in four otherwise identical models with differing initial conditions. These four models adopt the 
initial conditions shown in Figure [7] Upper left, upper right, lower left, and lower right panels show models GCM5, GCM6, GCM7, and 
GCM8,initialized from a rest state, an eastward decaying jet, a westward decaying jet, and an eastward barotropic jet, respectively. Unlike 
Figure [8] this figure shows snapshots at 1042 days, giving a sense of the instantaneous variability in the zonal-mean zonal wind. In a time 
average, this variability averages out and the different models would look even more similar. 
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Fig. 12. — Temperature (color scale, in K) and winds (arrows) for the same four models shown in Figure [TT] Panels show models GCM5 
(top left), GCM6 (top right), GCM7 (bottom left), and GCM8 (bottom right), initialized with a rest state, eastward decaying jet, westward 
decaying jet, and eastward barotropic jet, respectively. These are snapshots at 1042 days, and the slight differences between panels gives a 
sense of the time variability. In a time average, the different models look even more similar. 



ity in their models. By comparison, our models are set up 
so that the thermal forcing and horizontal entropy gradi- 
ents are weak at the lower boundary; this helps to avoid 



such lower-boundary instabilities, which are artificial in 
the context of a gas giant. 

Overall, our results indicate that specification of initial 
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Fig. 13. — Zonal-mean zonal wind at four different times for a single model, GCM8. Panels give snapshots at 833 days (top left), 903 
days (top right), 972 days (bottom left), and 1042 days (bottom right). The differences result from time variability. These differences are 
at least as large as the model-to-model differences shown in Figure [TT] indicating that those differences result from time variability rather 
than sensitivity of the statistical steady state to initial conditions. 



conditions is not a source of uncertainty in atmospheric 
circulation models of hot Jupiters, at least over the pa- 
rameter range explored here. This supports the contin- 
ued use of hot- Jupiter GCMs for understanding dynam- 
ical mechanisms, explaining observations, and making 
predictions to help guide future observations. That said, 
our results, as well as those of numerous previous pub- 
lications, show that details of the radiative forcing and 
frictional damping significantly affect the flow structure, 
including the qualitative dynamical regime, wind speeds, 
day- night temperature differences, and longitudinal off- 
sets of any hot or cold regions. State-of-the-art GCMs 
no w exist that include de tailed non-grey radiative trans- 
fer ( Showman et a l. 2009) as well as simpler, faster, gray 
treatments fH eng et al. 2 011a| |Rauscher fc Menou|2012 
'Perna et al. ^2012j) . By comparison, our understanding o: 



how to specify frictional drag is less well developed, and 



areas such as inclusion of clouds, sub-gridscale parame- 
terizations of turbulent mixing, and coupling to chem- 
istry have received little attention. Continued model de- 
velopment in these areas, and comparison of such mod- 
els to observations, should improve our ability to discern 
the physical and dynamical regimes of these fascinating 
worlds. 



We thank the International Summer Institute for Mod- 
eling in Astrophysics (ISIMA), and its organizers, par- 
ticularly Pascale Garaud and Doug Lin, for a stimulat- 
ing environment in which to begin this project. We also 
thank Yuan Lian, Lorenzo Polvani, and Ray Pierrehum- 
bert for useful discussions. 



REFERENCES 



Adcroft, A., Campin, J.-M., Hill, C, & Marshall, J. 2004, 

Monthly Weather Review, 132, 2845 
Arakawa, A., Sz Lamb, V. 1977, Methods in Computational 

Physics, 17, 173 
Budyko, M. I. 1969, Tellus, 21, 611 

Charbonneau, D., Knutson, H. A., Barman, T., Allen, L. E., 
Mayor, M., Megeath, S. T., Queloz, D., & Udry, S. 2008, ApJ, 
686, 1341 

Cho, J. Y.-K., & Polvani, L. M. 1996, Science, 8, 1 
Cooper, C. S., & Showman, A. P. 2005, ApJ, 629, L45 
— . 2006, ApJ, 649, 1048 

Cowan, N. B., Agol, E., & Charbonneau, D. 2007, MNRAS, 379, 
641 

Dobbs-Dixon, I., Cumming, A., & Lin, D. N. C. 2010, ApJ, 710, 
1395 

Dobbs-Dixon, I., & Lin, D. N. C. 2008, ApJ, 673, 513 
Dowling, T. E., Sz Ingersoll, A. P. 1989, Journal of Atmospheric 
Sciences, 46, 3256 



Hack, J. J., Sz Jakob, R. 1992, Description of a global shallow 

water model based on the spectral transform method, Tech. 

rep., National Center for Atmospheric Research Technical note 

NCAR/TN-343+STR, Boulder, CO 
Harrington, J., Hansen, B. M., Luszcz, S. H., Seager, S., Deming, 

D., Menou, K., Cho, J. Y.-K., & Richardson, L. J. 2006, 

Science, 314, 623 
Harrington, J., Luszcz, S., Seager, S., Deming, D., Sz Richardson, 

L. J. 2007, Nature, 447, 691 
Held, I. M., Sz Suarez, M. J. 1994, Bulletin of the American 

Meteorological Society, vol. 75, Issue 10, pp.1825-1830, 75, 1825 
Heng, K., Frierson, D. M. W., Sz Phillipps, P. J. 2011a, MNRAS, 

418, 2669 

Heng, K., Menou, K., & Phillipps, P. J. 2011b, MNRAS, 413, 2380 
Iro, N., Bezard, B., Sl Guillot, T. 2005, A&A, 436, 719 
Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., Sz 

Megeath, S. T. 2008, ApJ, 673, 526 
Knutson, H. A., et al. 2007, Nature, 447, 183 



13 



Lewis, N. K., Showman, A. P., Fortney, J. J., Marley, M. S., 

Freedman, R. S., & Lodders, K. 2010, ApJ, 720, 344 
Li, J., Sz Goodman, J. 2010, ArXiv e-prints 
Marley, M. S., & McKay, C. P. 1999, Icarus, 138, 268 
Mayor, M., Sz Queloz, D. 1995, Nature, 378, 355 
Menou, K. 2012, ApJ, 745, 138 
Menou, K., & Rauscher, E. 2009, ApJ, 700, 887 
Perna, R., Heng, K., & Pont, F. 2012, ApJ, 751, 59 
Perna, R., Menou, K., & Rauscher, E. 2010, ApJ, 719, 1421 
Pierrehumbert, R. T. 2010, Principles of Planetary Climate 

(Cambridge University Press) 
Polvani, L. M., Waugh, D. W., Sz Plumb, R. A. 1995, Journal of 

Atmospheric Sciences, 52, 1288 
Rauscher, E., Sz Menou, K. 2010, ApJ, 714, 1334 
— . 2012, ApJ, 745, 78 

Scott, R. K., & Polvani, L. 2007, J. Atmos. Sci, 64, 3158 
Scott, R. K., Sz Polvani, L. M. 2008, Geophys. Res. Lett., 35, 
L24202 

Showman, A. P. 2007, J. Atmos. Sci., 64, 3132 
Showman, A. P., Cooper, C. S., Fortney, J. J., & Marley, M. S. 
2008, ApJ, 682, 559 



Showman, A. P., Fortney, J. J., Sz Lewis, N. 2012, ApJ, 

submitted to ApJ 
Showman, A. P., Fortney, J. J., Lian, Y., Marley, M. S., 

Freedman, R. S., Knutson, H. A., Sz Charbonneau, D. 2009, 

ApJ, 699, 564 
Showman, A. P., & Guillot, T. 2002, A&A, 385, 166 
Showman, A. P., & Polvani, L. M. 2010, Geophys. Res. Lett., 37, 

18811 

— . 2011, ApJ, 738, 71 

Thrastarson, H. T., Sl Cho, J. 2010, ApJ, 716, 144 

Thrastarson, H. T., & Cho, J. Y. 2011, ApJ, 729, 117 

Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics: 

Fundamentals and Large-Scale Circulation (Cambridge Univ. 

Press, Cambridge, UK) 
Watkins, C, & Cho, J. 2010, ApJ 
Wright, J. T., et al. 2011, PASP, 123, 412 



